      PROGRAM QMC_TB
!========+=========+=========+=========+=========+=========+=========+=
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ 
! ----------  ----------  ----------  ----------  ----------  ----------
!     QMC + TB  (tight-binding)
!
!  Comments:
!  The initial Weiss GF is in file  "g0init".
!  Partial DOSes should be in pdos.in file.
!  Nd - Number of Degenerate spin-orbitals (Nd = 2N only diagonal!!!)
!  Nf - Number of auxilary Ising Fields (Nf = Nd(Nd-1)/2)
!
!  Last modified @ June 27 2k3 by Viktor Oudovenko   
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ 
      include 'param.dat'

! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI      
      include 'mpif.h'
#endif      
! ----------  ----------  ----------  ----------  ----------  ----------

      INTEGER  MyID,ierr,numprocs,MASTER,ibeta, fh_info, NDOS

      COMPLEX*16 G(nlm,nom,ns),     G0(nlm,nom,ns),
     &           Sig(nlm,nom,ns),   Sig2(nlm,nom,ns),
     &           delta(nlm,nom,ns), Gom(nom)

      COMPLEX*16 omes, ctmp, ctmp2, ctmp4, ctmpt, omes0, cunity
     
      REAl*8       Gt(L), GSUN(L)
     &            ,Gt0(L,Nd),   Gtmp(L,Nd),   Gmp(L,Nd)
     &            ,Gstat(L,Nd), Gstatm(L,Nd), Gav(L,Nd)
     &            ,zef(nlm,2),  U(Nf), Umat(Nd,Nd)
     &            ,omega(nom), omr(nom)
     &            ,nn(Nf), nnmat(Nd,Nd), nnt(Nf), vn(L,Nf)
     &            ,dos(nlm,nom,ns)
     &            ,FMA(Nd), SMA(Nd)
     &            ,FDA1(Nd), FDA2(Nd) 
     &            ,FM_DOS(nlm,ns), SM_DOS(nlm,ns), Vk2(nlm,ns)
                 
      REAl*8       occ, DelGF, DelGF0, X,Y, dw, dwold

      CHARACTER*1  IMI, CBETA1
      CHARACTER*2  CBETA
      CHARACTER*4  CU
      CHARACTER*6  CMU
      CHARACTER*19 FILEN

      LOGICAL TB, SC_TEST, SUN, UINFO, LSC, FD_SIMPLE

! ----------  ----------  ----------  ----------  ----------  ----------
      COMMON/varmpi/  MyID,ierr,numprocs,MASTER
      COMMON/param/   domega,dtime,omega,temp
      COMMON/ising/   vn
! ----------  ----------  ----------  ----------  ----------  ----------
      MyID = 0
      numprocs = 1

! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI     
! getting in contact with MPI
      CALL MPI_INIT( ierr )
      CALL MPI_COMM_RANK( MPI_COMM_WORLD, MyID, ierr )
      CALL MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
#endif     
! ----------  ----------  ----------  ----------  ----------  ----------

      MASTER = 0
	Ncpu = numprocs
!------------> Global constant:
      Zero = 0.d0
      One  = 1.d0
      Two  = 2.d0

      IF (IWMAX.GE.2**16) PRINT*, 'STOP : IWMAX IS TOO LARGE'
 
      UINFO = .FALSE.  
! ----------  ----------  ----------  ----------  ----------  ----------
	IF (MyID .EQ. MASTER) THEN  ! MASTER

	  IF (.NOT.UINFO)	THEN
      	fh_info = 500
          OPEN(fh_info,FILE = 'usr_info',status = 'unknown')
          CALL userinfo(fh_info)
	  ENDIF  !(.NOT.UINFO)	


      OPEN(1, file = 'inp', STATUS = 'UNKNOWN')

      READ(1,*) uu,    uj
      READ(1,*) Beta,  Nsweep
      READ(1,*) ef0,   def, znm
      READ(1,*) istart,    nscf,  iout
      READ(1,*) alpha, small
      READ(1,*) NewFour, FD_SIMPLE
      READ(1,*) TB, LSC

!------------> INPUT EXPLANATION --------------------------------------
!     uu,uj  - Hubbard parameters: U and J
!
!     Beta   - inverse temperature
!     Nsweep - # of Mote-Carlo sweeps
!
!     ef0    - the chemical potential
!     def    - Energy step for Fermy Energy search 
!            = 0 the chem potential is fixed and program finds doping
!           != 0 the doping is fixed (znum) and program searches for ef
!     znum   - Guess for doping 
!
!     istart = 0 (Start:not READ Sigma) 
!            = 1 (READ Sig but not WRITE) 
!            = 2 (READ and re-WRITE Sig) 
!     nscf   - max. # of iterations
!     iout   - control of output files
!            - 0 keep only one last iteration in output files
!
!     alpha  - small mixing coefficient (of the order of 0.05-0.5)
!     small  - Selfconsistency threshold for Delta(Sig)/Average(Sig) 
!
!     NewFour- # of Fourier Transformation (FT) routine used
!            = 1 the best scheme (the 1-st derivatives to be provoded)     
!            = 2 a good scheme (the moments to be provoded)
!            = 3 worst scheme (natural spline, nothing is known)
!     FD_SIMP- in scheme N1 of FT provides the first derivatatives 
!              using extrapolation (LOGICAL)
!
!     TB     = F ->  Semicircular DOS  (LOGICAL)
!            = T ->  TB DOS in file pdos.in (must be provided) (LOGICAL)
!     LSC    - reserved for external k-summations (LOGICAL)
! ----------  ----------  ----------  ----------  ----------  ----------

      temp = One/Beta
      dtau = Beta/dfloat(L)
      t = 0.5D0   !  
      t2 = t**2     
	Vk = t**2 !it should be calculated from: int -1/pi G0 dw = \pi Vk  

      IF (UJ.Eq.0) SUN = .TRUE.

      PRINT*,'QMC: multi-band Hubbard model N = ',nlm
      PRINT*,'Number of CPUs  = ', numprocs
      PRINT*,'USERINFO, SU(N) = ', UINFO, SUN
      PRINT*,' '
!
      CALL matsub
      CALL vertex(uu,uj,u,umat)
!
      PRINT*,'------------------------------------------------------- '
      PRINT*,'U,  J       = ', real(uu),real(uj)
      PRINT*,'U(ij)       = '
      WRITE(6,*)(U(IF),IF = 1,Nf)
      PRINT*,' '
      PRINT*,'Input  Ef0  = ', ef0
      PRINT*,' '
      PRINT*,'Beta, dtau  = ',real(beta),real(dtau)
      PRINT*,'T(eV,K)     = ',real(temp),real(temp*11605)
      PRINT*,'dT,   dOm   = ',real(dtime),real(domega)
      PRINT*,'             '
      PRINT*,'istart      = ', istart
      PRINT*,'New FOURIER = ',NewFOUR
      PRINT*,'TB          = ',TB
      PRINT*,'t           = ',t
      PRINT*,'Starting number of electrons ZNM = ',znm
	PRINT*,'------------------------------------------------------- '
  
! ----------  ----------  ----------  ----------  ----------  ----------

! Preparation of output filenames which should contain information about
! interaction, U, temperature, beta, the chemical potential, ef.

      IF (UU.ge.10) then
	  WRITE(CU ,3209)  UU
	ELSE
	  WRITE(CU ,3210)  UU
	ENDIF
      
      ibeta = beta
	IF (beta.ge.10) then
              WRITE(CBETA,3213)  ibeta
       ELSE
              WRITE(CBETA1,3211)  ibeta
                    CBETA = "0"//CBETA1
       ENDIF

      IF0u = ef0      
      IF (abs(ef0).ge.10) then
	   WRITE(CMU,3209)  abs(ef0)    
	ELSE
	   WRITE(CMU,3210)  abs(ef0)    
	ENDIF
	
      IF (IF0u.le.0) THEN
         FILEN = "_"//"b"//CBETA//"_"//CU//"_"//"m"//CMU
      ELSE
         FILEN = "_"//"b"//CBETA//"_"//CU//"_"//"p"//CMU
      ENDIF

3209  FORMAT(F5.2)
3210  FORMAT(F4.2)
3212  FORMAT(F5.3)
3211  FORMAT(I1)
3213  FORMAT(I2)

!     Openning necessary files
      OPEN(2, file = 'Gtau'//FILEN,STATUS = 'UNKNOWN')
	IF (istart .EQ. 0)  THEN
      OPEN(10,file='Sig'//FILEN, STATUS = 'UNKNOWN',form='unFORMATted')
	ELSE
      OPEN(10,file='Sig',        STATUS = 'UNKNOWN',form='unFORMATted')
	ENDIF
! ----------  ----------  ----------  ----------  ----------  ----------
      OPEN(13,file = 'g0init'//FILEN,  STATUS = 'UNKNOWN')
      OPEN(14,file = 'g0init',		 STATUS = 'UNKNOWN')
      OPEN(22,file = 'G0tau'//FILEN,   STATUS = 'UNKNOWN')
      OPEN(44,file = 'greenw'//FILEN,  STATUS = 'UNKNOWN')
      OPEN(46,file = 'greenw0'//FILEN, STATUS = 'UNKNOWN')
      OPEN(54,file = 'sigma'//FILEN,   STATUS = 'UNKNOWN')
      OPEN(67,file = 'moments'//FILEN, STATUS = 'UNKNOWN')
      CLOSE(1)
! ----------  ----------  ----------  ----------  ----------  ----------

      IF (TB)  THEN
!      IF (TB .AND. .NOT.LSC )  THEN

      OPEN(55,file = 'pdos.in',   STATUS = 'UNKNOWN')
      OPEN(56,file = 'pdos.check',STATUS = 'UNKNOWN')
 
!-----------  Read LDA partial DOSes --------------- 
      IF ( SUN ) THEN

        DO i = 1,nom 
          READ(55,*,END = 55)ener,doss 
          omr(I) = ener 
 
          IF (i .EQ. 1) emin = ener 
	      DO im =1, nlm
              dos(im,i,1) = doss 
	      ENDDO
        ENDDO 

	ELSE ! ( SUN )

        DO i = 1,nom 
          READ(55,*,END = 55)ener, ((dos(im,i,is), im = 1,nlm),is=1,ns)
          omr(I) = ener 
          IF (i .EQ. 1) emin = ener 
        ENDDO 

	ENDIF !( SUN )
55    continue 
      emax = ener 
      ndos = i-1 
      dw = (emax-emin)/float(ndos-1) 
 
      WRITE(6,*) 
      WRITE(6,"(' DOS: emin,emax,ndos',2(3x,f8.4),3x,i5)")emin,emax,ndos 
      WRITE(6,*) 
      WRITE(6,"(' QMC: Nsweep*Numprocs',2x,i7,'* ',i3)") nsweep,numprocs 
      WRITE(6,*) 
 
 
!-----------  Checking frequency mesh --------------- 
!      DO I = 2, ndos-1 
!         IF ( (omr(i)-omr(i-1)).NE.(omr(i+1)-omr(i)) )  
!     &   STOP 'Something wrong with frequency mesh, plz, check!' 
!      ENDDO 
 
!-----------  Calculation of the first DOS moment ---- 

      PRINT *, 'The first moment for DOS = '
	DO is = 1,ns
	  DO im = 1,nlm
          sum = 0
	    sum2= 0
          DO I = 1, ndos           
            sum = sum  + dos(im,i,is)*dw*omr(i)
	      sum2= sum2 + dos(im,i,is)*dw*omr(i)**2
            WRITE(56,*) omr(i),dos(im,i,is)
          ENDDO 
          FM_DOS(im,is) = sum
	    SM_DOS(im,is) = sum2
	    Vk2   (im,is) = sum2-sum**2
	    PRINT *, is, im, sum, sum2, Vk2(im,is)
        ENDDO
	ENDDO    

      CLOSE(56)

	ELSE !(TB)

	DO is = 1,ns
	  DO im = 1,nlm
          FM_DOS(im,is) = 0.0
        ENDDO
	ENDDO    

	ENDIF !(TB)

	GOTO 113

      DO i = 1,nom 
         omr(I) = omr(I) - sum
      ENDDO 


      PRINT *, 'The first moment for DOS = '
	DO is = 1,ns
	  DO im = 1,nlm
          sum = 0
	    sum2= 0
          DO I = 1, ndos           
            sum = sum  + dos(im,i,is)*dw*omr(i)
	      sum2= sum2 + dos(im,i,is)*dw*omr(i)**2
            WRITE(56,*) omr(i),dos(im,i,is)
          ENDDO 
          FM_DOS(im,is) = sum
	    SM_DOS(im,is) = sum2
	    Vk2   (im,is) = sum2-sum**2
	    PRINT *, is, im, sum, sum2, Vk2(im,is)
        ENDDO
	ENDDO    

113   CONTINUE
! ----------  ----------  ----------  ----------  ----------  ----------
!------------> Start
      ef = Ef0 
      IF (istart .EQ. 0) then
        DO 5 is = 1,ns
        DO 5 io = 1,nom
        DO 5 im = 1,nlm
 5      Sig(im,io,is) = Dcmplx(0.0d0,0.0d0)
      ELSE  ! istart
        REWIND 10
        READ(10)(((Sig(im,io,is),im = 1,nlm),io = 1,nom),is = 1,ns)
      ENDIF ! istart
	IF ( istart.ne.0 )  THEN
	CLOSE(10)
      OPEN(10,file='Sig'//FILEN,STATUS='UNKNOWN',form='unFORMATted')
	ENDIF

!------------> G0(tau)
      REWIND(14)
      DO I = 1,L
        READ(14,*)(gt0(I,m),m = 1,nlm)
      ENDDO

      DO I = 1,L
        WRITE(13,'(I4,$)') I
        WRITE(13,'(f22.16,$)') (Gt0(i,m),m = 1,nlm)
        WRITE(13,*)
      ENDDO
        WRITE(13,'(I4,$)' )L+1
        WRITE(13,'(f22.16,$)')((-Gt0(1,m)-1.d0),m = 1,nlm)
      
	CALL flush(13)
! ----------  ----------  ----------  ----------  ----------  ----------
      ENDIF  ! MASTER
! SEND the data we have READ from process 0 (MASTER) to the other processes.
      Numf = Nf
      Numg = Nd*L

! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI   

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_BCAST(u, Numf, MPI_DOUBLE_PRECISION, MASTER,
     &                                            MPI_COMM_WORLD, ierr)

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_BCAST(small, 1, MPI_DOUBLE_PRECISION, MASTER,
     &                                            MPI_COMM_WORLD, ierr)

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_BCAST(dtau,  1, MPI_DOUBLE_PRECISION, MASTER,
     &                                            MPI_COMM_WORLD, ierr)

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_BCAST(nscf,  1, MPI_INTEGER, MASTER,
     &                                            MPI_COMM_WORLD, ierr)

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_BCAST(Nsweep,1, MPI_INTEGER, MASTER,
     &                                            MPI_COMM_WORLD, ierr)
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )

!------------> random number seed Idum for each procs.
      IF (MyID .EQ. MASTER) then  ! MASTER
 
      Idum = 135749
      Inum = 989769

        DO i = 1,numprocs-1
        Idum = int(ranw(Idum)*Inum) + 1
        CALL MPI_SEND(Idum, 1, MPI_INTEGER, i,99, MPI_COMM_WORLD, ierr)
        ENDDO
        Idum = int(ranw(Idum)*Inum) + 1
      ELSE
       CALL MPI_RECV(Idum, 1, MPI_INTEGER, MASTER,
     &           99, MPI_COMM_WORLD, status, ierr)
!
      endif  ! MASTER
#endif
! ----------  ----------  ----------  ----------  ----------  ----------      
! __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __
!
      DO im = 1,Nlm
       DO is = 1,Ns
        DO il = 1,L
         gt(il) = Gt0(il,im+(is-1)*Nlm)
        ENDDO

          IF     ( NewFOUR .EQ. 1 .OR. NewFOUR .EQ. 2) THEN
            CALL FDfit( Gt,dtau, FD1,FD2 )
		  CALL nfourier1( FD1,FD2,gt,gom )
	
          ELSE
            CALL nfourier(gt,gom)


          ENDIF ! NewFOUR

        DO io = 1,Nom
          G0(im,io,is) = gom(io)
        ENDDO
       ENDDO
      ENDDO

! __  __  __  __  __  __ Self-consistent run  __  __  __  __  __  __  __ 
      dS = 1.
      aS = 1.
      dSaS = dS/aS
! AL's choice  (not tested properly)
	Eft = 0       
	Xmu = Ef  
! The original one	    
	Eft = ef       
	Xmu = 0      
! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI   
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )                  
      CALL MPI_BCAST(Eft, 1, MPI_DOUBLE_PRECISION, master,      
     &               MPI_COMM_WORLD, ierr)                      
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )                  
      CALL MPI_BCAST(Xmu, 1, MPI_DOUBLE_PRECISION, master,      
     &               MPI_COMM_WORLD, ierr)                      
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )                  
#endif      
! ----------  ----------  ----------  ----------  ----------  ----------
      DO 1 ISCF = 1,NSCF  

      DO in = 1,Nf
          nn(in) = 0.d0
          nnt(in) = 0.d0
      ENDDO
! ----------  ----------  ----------  ----------  ----------  ----------
!------------> output
      IF (dSaS.le.small) goto 2   ! now it works
      IF (MyID .EQ. MASTER) then  ! MASTER
!
      PRINT*,'****** Iteration Number = ',iscf
      IF (iout .EQ. 0) REWIND 2
      IF (iout .EQ. 0) REWIND 3
      IF (iout .EQ. 0) REWIND 22
!      IF (iout .EQ. 0) REWIND 33

!------------>  READ in initial green's function 
      REWIND(13)
      DO i = 1,L
        READ(13,*)itmp,(gt0(i,m),m = 1,nlm)
      ENDDO
	      
      DO i = 1,L
        DO m = 1,nlm
           gt0(i,m+nlm) = gt0(i,m) 
        ENDDO
      ENDDO
!
       DO il = 1,L
        DO id = 1,Nd
         gtmp(il,id) = 0.d0
        ENDDO
       ENDDO
! ----------  ----------  ----------  ----------  ----------  ----------
      ENDIF  ! MASTER

! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI   
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_BCAST(gt0, Numg, MPI_DOUBLE_PRECISION, MASTER,
     &               MPI_COMM_WORLD, ierr)
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
#endif      
! ----------  ----------  ----------  ----------  ----------  ----------
!                          Q   M  C   
	 CALL qmc(Gt0,gmp,Xmu,dtau,u,iscf,idum,acc,nrat,Nsweep,nnt)
!
! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI  
!------------> Sum-up Gtmp from all processors 
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_REDUCE( gmp, gtmp, Numg, MPI_DOUBLE_PRECISION,
     &                MPI_SUM, MASTER, MPI_COMM_WORLD, ierr )

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_REDUCE( acc, acctot, 1, MPI_DOUBLE_PRECISION,
     &                MPI_SUM, MASTER, MPI_COMM_WORLD, ierr )

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_REDUCE( nrat, nratot, 1, MPI_INTEGER,
     &                MPI_SUM, MASTER, MPI_COMM_WORLD, ierr )

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_REDUCE( nnt, nn, Nf, MPI_DOUBLE_PRECISION,
     &                MPI_SUM, MASTER, MPI_COMM_WORLD, ierr )

      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
#endif      
! ----------  ----------  ----------  ----------  ----------  ----------
!
      IF ( MyID .EQ. MASTER ) then
      WRITE(*,*) 'acc.rate:',acctot/dfloat(numprocs)
      WRITE(*,*) 'acc.rate:',acc
      WRITE(*,*) 'Total number of neg det:',nratot
      WRITE(*,*) 'Total number of neg det:',nrat
      WRITE(*,*) '_____________________________________________________'
      WRITE(*,*) '--------------------- nn ----------------------------'

        IF (numprocs .EQ. 1) THEN
	        DO i = 1,L
       		  DO m = 1,Nd
       		   Gtmp(i,m) = Gmp(i,m)
       		  ENDDO
      		ENDDO

        	    acctot = acc
      		nratot = nrat

	      	DO m = 1,Nf
       		 nn(m) = nnt(m)
      		ENDDO
		    PRINT *, 'N CPUs = ', numprocs
        ELSE
	        DO i = 1,L
       		  DO m = 1,Nd
       		   Gtmp(i,m) = Gtmp(i,m)/float(numprocs)
       		  ENDDO
      		ENDDO
        ENDIF

! the double occupanies in SU(N) case
      docc = 0.D0
      DO m = 1,Nf
        WRITE(*,*) m, nn(m)
	docc = nn(m) +docc
      ENDDO 
	docc = docc/float(Nf)
      WRITE(*,*) 'docc = ', docc
      WRITE(*,*) '---------------------------------------------------'

! the double occupanies
      ij=0.0
      do i=1,Nd-1
        do j=i+1,Nd          
	    ij=ij+1
	    nnmat(i,j)=nn(ij)
        enddo
      enddo
      do i = 1,Nd
	  nnmat(i,i) = 0.0
        do j = i,Nd
          nnmat(j,i) = nnmat(i,j)
        enddo
      enddo

!------------> Nonmagnetic case - average Gtau
      IF (Ns .EQ. 1) then
        DO i = 1,L
          DO m = 1,Nlm
            Gtmp(i,m) = (Gtmp(i,m)+Gtmp(i,m+nlm))/TWO
            Gtmp(i,m+nlm) = Gtmp(i,m)
          ENDDO 
        ENDDO 
      ENDIF
!
! SU(N) GF
      IF ( SUN ) THEN
        DO i = 1,L
 	    GSUN(i) = 0.0
          DO m = 1,Nd
            GSUN(i) = Gtmp(i,m)+GSUN(i)
          ENDDO
          GSUN(i) = GSUN(i)/float(Nd)
          DO m = 1,Nd
            Gtmp(i,m) = GSUN(i)
            Gav (i,m) = GSUN(i)
            Gstatm(i,m) = 0.0
          ENDDO
        ENDDO
      ENDIF

      ENDIF  ! MASTER

! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI      
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
! ----------  ----------  ----------  ----------  ----------  ----------
! COMPUTATION OF Delta_g for MaxEnt program
! ----------  ----------  ----------  ----------  ----------  ----------
      IF (MyID .EQ. MASTER) then  ! MASTER
      Numg = Nd*L
        DO i = 1,numprocs-1
        CALL MPI_SEND(Gav,Numg,MPI_DOUBLE_PRECISION,i, 
     &                 99, MPI_COMM_WORLD, ierr)
        ENDDO
      ELSE
      Numg = Nd*L
      CALL MPI_RECV(Gav, Numg, MPI_DOUBLE_PRECISION, MASTER,
     &           99, MPI_COMM_WORLD, status, ierr)
!
      endif  ! MASTER
! __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __ 
        DO i = 1,L
         DO m = 1,Nd
          Gstat(i,m) = ( Gmp(i,m)-Gav(i,m) )**2
         ENDDO
        ENDDO
! ----------  ----------  ----------  ----------  ----------  ----------
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_REDUCE( Gstat, Gstatm, Numg, MPI_DOUBLE_PRECISION,
     &                MPI_SUM, MASTER, MPI_COMM_WORLD, ierr )
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
! ----------  ----------  ----------  ----------  ----------  ----------
#endif  
! ----------  ----------  ----------  ----------  ----------  ----------

      IF ( MyID .EQ. MASTER ) then

! Statistics for maxent program
	IF (NCPU.GT.1) THEN
        DO i = 1,L
         DO m = 1,Nd
           Gstatm(i,m) = DSQRT(Gstatm(i,m)/float((ncpu-1)*ncpu))
         ENDDO
        ENDDO
	ELSE
	DO i = 1,L
         DO m = 1,Nd
           Gstatm(i,m) = 0.0
         ENDDO
      ENDDO
	ENDIF  !(NCPU.GT.1)
! ----------  ----------  ----------  ----------  ----------  ----------
!------------> G(tau) OUTPUT 

	IF ( ns .EQ. 1) THEN

      DO i = 1,L
       WRITE(2,'(F16.10,$)')dtau*(i-1),((Gtmp(i,m),Gstatm(i,m)),m=1,nlm)
       WRITE(2,*)
      ENDDO
       WRITE(2,'(F16.10,$)') beta,
     & (((-Gtmp(1,m)-1.d0),Gstatm(1,m)),m = 1,nlm)
       WRITE(2,*)


	ELSE  ! ( ns .EQ. 1) 

      DO i = 1,L
       WRITE(2,'(F16.10,$)')dtau*(i-1),
     & (((Gtmp(i,m+nlm*(is-1)),Gstatm(i,m+nlm*(is-1))),m=1,nlm),is=1,ns)
       WRITE(2,*)
      ENDDO
       WRITE(2,'(F16.10,$)') beta,
     &((((-Gtmp(1,m+nlm*(is-1))-1.d0),
     &   Gstatm(1,m+nlm*(is-1))),m=1,nlm),is=1,ns)
       WRITE(2,*)


	ENDIF ! ( ns .EQ. 1)
! ----------  ----------  ----------  ----------  ----------  ----------
!------------> Fourie to Omega  (Direct)
      IF ( SUN ) THEN
       IF ( iscf .EQ. 1 ) THEN
          occ = znm/Two/DFLOAT(Nlm)  
	    docc = 0
       ELSE
            occ = zf/ns/dFLOAT(Nlm)
       ENDIF
!	
          IF      ( NewFOUR .EQ. 1 ) THEN
	      CALL FDeSSUN( omega,eft,G0,G,FM_DOS,UU,occ,docc, FDA1,FDA2)

	    ELSE	    
	      CALL moment(Nlm,UU,ef,occ,docc,FM_DOS(1,1),Vk2(1,1),FM,SM)
            WRITE(6,*)'Moments 2,    3,      occ,        docc '
            WRITE(6,'(4F10.5)') FM, -SM, occ, docc
	    
		ENDIF


      ELSE   ! ( SUN ) 

          IF      ( NewFOUR .EQ. 1 ) THEN
		  IF ( iscf .EQ. 1)  THEN
  	        DO m = 1,Nlm
                zef(m,1) = znm/Two/DFLOAT(Nlm)  
	          zef(m,2) = zef(m,1)   
              ENDDO	  
	      ENDIF 
            CALL FDeS( omega,eft,G0,G,FM_DOS,Umat,zef,nnmat, FDA1,FDA2)

	    ELSE	    
	      CALL moments(Umat, ef, zef, nnmat, FM_DOS, Vk2, FMA, SMA)

		ENDIF

      ENDIF  ! ( SUN ) 
! ----------  ----------  ----------  ----------  ----------  ----------
       DO im = 1,Nlm
        DO is = 1,Ns
		 ims = im + (is-1)*Nlm
         DO il = 1,L
           gt(il) = Gtmp(il,im+(is-1)*Nlm)
         ENDDO

!      FT scheme N1 (using the first derivatives)
          IF      ( NewFOUR .EQ. 1) THEN
	      FD1 = FDA1(ims)
	      FD2 = FDA2(ims)
            IF ( FD_SIMPLE ) CALL FDfit(Gt,dtau, FD1,FD2)
		  CALL nfourier1( FD1,FD2,gt,gom )

!      FT scheme N2 (using moments)
	    ELSE IF ( NewFOUR .EQ. 2) THEN
		  IF (.NOT. SUN) THEN
	        FM = FMA(ims)
	        SM = SMA(ims)
	      ENDIF
            CALL nfourier2(One, FM, -SM, gt,gom)


          ELSE
!      FT scheme N3 (using the natural spline, old FT)
            CALL nfourier(gt,gom)
            PRINT *, 'OLD FOURIER '

          ENDIF

         DO io = 1,Nom
           G(im,io,is) = gom(io)
         ENDDO
        ENDDO
       ENDDO
!
! Output GF(w)
      IF (mod(iscf,6) .EQ. 0) REWIND 44   
      DO io = 1,nom
       WRITE(44,'(F16.9,$)')omega(io),(G (im,io,1),im = 1,nlm)
       WRITE(44,*)
      ENDDO
	CALL flush(44)
! ----------  ----------  ----------  ----------  ----------  ----------


      IF (TB)  THEN
!------------> DMFT Self-energy:  Sigma
      DO im = 1,Nlm
       DO is = 1,Ns
        DO io = 1,Nom
           Sig2(im,io,is) = One/G0(im,io,is)-One/G(im,io,is)
        ENDDO
       ENDDO
      ENDDO
!------------> mixing  Sig
	IF ( ISCF.GT.1.or.istart.ne.0 ) THEN
         CALL mix(Sig,Sig2,alpha,dS,aS)
	ELSE
         DO im = 1,Nlm
          DO is = 1,Ns
           DO io = 1,Nom
             Sig(im,io,is) = Sig2(im,io,is)
           ENDDO
          ENDDO
         ENDDO
	dS = 1
	aS = 1
	ENDIF
	dSTB = dS
	aSTB = aS


      DO 818 im = 1,nlm
      DO 818 is = 1,ns

!       "cunity" is necessary to ensure 1/omega limit when omega -> infinity
        cunity = 0.0
        DO i = 1,Ndos
           cunity = cunity+DOS(im,i,is)
        ENDDO
        cunity = cunity*dw
        PRINT *, '1 = ', cunity, im

      DO 818 io = 1,nom
         ctmp = 0.d0
         omes0 = ci*omega(io)+eft-Sig(im,io,is)
!
         DO i = 1,Ndos
            omes = omes0-omr(i)
            ctmp = ctmp+DOS(im,i,is)/omes
         ENDDO
!
         ctmp = ctmp / cunity 

         G(im,io,is) = ctmp*dw

         G0(im,io,is) = ( G0(im,io,is)+
     &   ONE/( ONE/G(im,io,is)+Sig(im,io,is) ) )/TWO

818    CONTINUE

! ---------- 
      ELSE ! (TB)
! ----------  
!      Starting to calculate G0 for semicircular dos
      PRINT*, 'Calculation Semicircular DOS', iscf

      DO 819 io = 1,nom
      DO 819 is = 1,ns
      DO 819 im = 1,nlm
         G0(im,io,is) = ( G0(im,io,is)+  
     &   One/(ci*omega(io)+eft-G(im,io,is)*t2 ) )/TWO

         Delta(im,io,is) = G(im,io,is)*t2  ! t^2*G (Hubbard)
!
819   CONTINUE

      ENDIF  !(TB)
! ----------  ----------  ----------  ----------  ----------  ----------

!  Output G0(w)
      REWIND 46   
      DO io = 1,nom
        WRITE(46,'(F16.10,$)') omega(io),(G0(im,io,1),im = 1,nlm)
        WRITE(46,*)
      ENDDO
!------------> Fourier to Time

      DO im = 1,Nlm
        DO is2 = 1,2
        
	    is = is2
          IF (Ns .EQ. 1) is = 1
            efr = -( eft )
	  
	      DO io = 1,Nom
              gom(io) = G0(im,io,is)
            ENDDO

            IF ( NewFour .EQ. 1  .OR. NewFour .EQ. 2) THEN
              CALL invfourierNEW(gom,gt,efR)
            ELSE
              CALL invfourier(gom,gt)
            ENDIF

         DO il = 1,L
          gt0(il,im+(is2-1)*Nlm) = gt(il)
        ENDDO

       ENDDO  
      ENDDO

! ----------  ----------  ----------  ----------  ----------  ----------
      REWIND(13)
      DO i = 1,L
        WRITE(13,'(I4,$)'    ) I
        WRITE(13,'(F16.11,$)') (Gt0(i,m),m = 1,nlm)
        WRITE(13,*)
      ENDDO
      WRITE(13,'(I4,$)'    ) L+1
      WRITE(13,'(F16.11,$)') ((-Gt0(1,m)-1.d0),m = 1,nlm)
      WRITE(13,*)
      
!------------> G0(tau)
      DO i = 1,L
        WRITE(22,'(F16.10,$)') dtau*(i-1),(Gt0(i,m),m = 1,nlm)
        WRITE(22,*)
      ENDDO
      WRITE(22,'(F16.10,$)') Beta,((-Gt0(1,m)-1.d0),m = 1,nlm)
      WRITE(22,*)
! ----------  ----------  ----------  ----------  ----------  ----------
!------------> Fermi  search
!------------> D-occupation from G(tau)
      DO m = 1,Nlm
        zef(m,1) = 1.d0+Gtmp(1,m)
	  zef(m,2) = zef(m,1)   !check
        IF (ns.gt.1) then
          zef(m,2) = 1.d0+Gtmp(1,m+Nlm)
        ENDIF
      ENDDO
         zfu = 0.d0
         zfd = 0.d0
      DO m = 1,nlm
         zfu = zfu+zef(m,1)
         IF (ns.gt.1) then
            zfd = zfd+zef(m,2)
         ENDIF
      ENDDO
      IF (ns.gt.1) then
         zf = zfu+zfd
      ELSE
         zf = zfu
      ENDIF

      PRINT*,'D-num-ef = ',zf
      PRINT*,'D-num-up  '
      WRITE(6,'(F10.5,$)') (zef(il,1),il = 1,nlm)
      WRITE(6,*)

      IF (ns.gt.1) then
          PRINT*,'D-num-dn  '
          WRITE(6,'(F10.5,$)') (zef(il,2),il = 1,nlm)
          WRITE(6,*)
          PRINT*,'num-ef,  Mag   total = '
          WRITE(6,'(2F10.5)') zf,zfu-zfd
      ELSE
          PRINT*,'num-ef, total = '
          WRITE(6,'(F10.5)') zf
      ENDIF
!------------> Ajusting Ef
      IF (NS .EQ. 1) THEN
        efn = ef+(znm-zf*2)*def
      ELSE
        efn = ef+(znm-zf)*def
      ENDIF
      PRINT*,'ef, ef-new = ',REAL(ef), REAL(efn)
!     ef = (1.d0-alpha)*ef + alpha*efn
      ef = efn
! ----------  ----------  ----------  ----------  ----------  ----------
	IF ( .NOT.TB )  THEN
!------------> DMFT Self-energy:  Sigma
       DO im = 1,Nlm
         DO io = 1,Nom
           DO is = 1,Ns
             Sig2(im,io,is) = One/G0(im,io,is)-One/G(im,io,is)
           ENDDO
         ENDDO
       ENDDO
!------------> mix  Sig
	IF ( ISCF.GT.1.or.istart.ne.0 ) THEN
       CALL mix(Sig,Sig2,alpha,dS,aS)
	ELSE
      DO im = 1,Nlm
        DO is = 1,Ns
          DO io = 1,Nom
             Sig(im,io,is) = Sig2(im,io,is)
          ENDDO
        ENDDO
      ENDDO
	ds = 1
	as = 1
	ENDIF

	ENDIF !(.NOT.TB)
! ----------  ----------  ----------  ----------  ----------  ----------
!------------> END SCF-loop
!          IF (iscf.LE.(nscf- 3).AND.IWMAX.GE.2**15) REWIND 54
!          IF (iscf.LE.(nscf- 5).AND.IWMAX.GE.2**14) REWIND 54
!          IF (iscf.LE.(nscf-20).AND.IWMAX.GE.2**13) REWIND 54
	    
		REWIND 54

          WRITE(54,*) '#iter = ',  iscf, ef
          DO io = 1,nom 
            WRITE(54,'(F16.9,$)') 
     &                           omega(io),(Sig(im,io,1),im = 1,nlm)
            WRITE(54,*)
          ENDDO

          PRINT*,'iscf,   ef,       dS/aS'
          PRINT*,'---------------------------------------------------'
	IF (aS .EQ. 0.d0) aS = 1.d0
        WRITE(6,1000) iscf, ef, dS/aS
1000  FORMAT(i4,8F10.4  )
      dSaS = dS/aS
	IF (TB)  dSaS = dSTB/aSTB
      ENDIF  ! MASTER

! ---------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI      
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_BCAST(dSaS, 1, MPI_DOUBLE_PRECISION, MASTER,
     &               MPI_COMM_WORLD, ierr)
#endif  
! ---------  ----------  ----------  ----------  ----------  ----------
1       CONTINUE
! __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __  __ 
2       CONTINUE

      IF ( MyID .EQ. MASTER ) then 

          IF ( istart.ne.1 ) then
            REWIND 10
            WRITE(10)(((Sig(im,io,is),im = 1,nlm),io = 1,nom),is = 1,ns)
          ENDIF
! ---------  ----------  ----------  ----------  ----------  ----------
      WRITE(67,*) 'beta = ', beta
      docc = 0.D0
      DO m = 1,Nf
        WRITE(67,*) m,nn(m)
        docc = nn(m) + docc
      ENDDO 
      docc = docc/float(Nf)
      WRITE(67,*) '--------------------------------------------------'
      WRITE(67,*) 'docc = ' , docc
      WRITE(67,*) (m, nn(m),m = 1,Nf)
      WRITE(67,*) '--------------------------------------------------'
! ----------  ----------  ----------  ----------  ----------  ----------

	IF ( .NOT.UINFO )	THEN
           CALL exitinfo(fh_info)
           CALL userinfo(fh_info)
           CLOSE(500)
	ENDIF ! (.NOT.UINFO)

      ENDIF  ! MASTER

! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI    
!------------> End MPI:
      CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
      CALL MPI_FINALIZE( ierr )
#endif  
! ----------  ----------  ----------  ----------  ----------  ----------

      END
! ----------  ----------  ----------  ----------  ----------  ----------
